Eulerian Statistically Preserved Structures in Passive Scalar Advection 



m 
o 
o 



in 

Q 

U 
d 



> 

o 
o 
m 
o 
m 
o 



> 

X 
S3 



Yoram Cohen 1 ' 2 Anna Pomyalov 1 and Itamar Procaccia 1 ' 2 
1 Dept. of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel. 
2 Dept. of Physics, the Chinese University of Hong Kong, Shatin, Hong Kong. 

We analyze numerically the time-dependent linear operators that govern the dynamics of Eulerian 
correlation functions of a decaying passive scalar advected by a stationary, forced 2-dimensional 
Navier-Stokes turbulence. We show how to naturally discuss the dynamics in terms of effective 
compact operators that display Eulerian Statistically Preserved Structures which determine the 
anomalous scaling of the correlation functions. In passing we point out a bonus of the present 
approach, in providing analytic predictions for the time-dependent correlation functions in decaying 
turbulent transport. 



I. INTRODUCTION 

The aim of this paper is to discuss the statistical 
physics of turbulent advection of passive scalars We 
are interested in scalar fields 8(r, t) which are advected 
by a velocity field u(r,t) such that together they solve 
the set of equations 



9 U , „v 



---Vp + vAu + f, (1) 
kAO + f . (2) 



In these equations p(r, t), v and re are the pressure field, 
the kinematic viscosity and the scalar diffusivity respec- 
tively. In this paper the Navier-Stokes equations will 
be always forced with a time-stationary forcing f(r,t). 
The scalar field will be forced or unforced (decaying), de- 
pending on our interests below. We focus on the case of 
high Reynolds number Re and high Peclet number Pe, 
i.e. v, re — ► 0, where the turbulence of the velocity field is 
fully developed and where both fields display a significant 
range of scaling behavior at scales that are sufficiently far 
from the forcing scales. 

A major theoretical question that had been answered 
recently has to do with the scaling properties of the 
correlation functions of the advected field 0, 13 De- 
fine the simultaneous many-point correlation function 
F^ N '(ri, . . . r N ) in the forced case by 

F^\r ll ...r N ) = (6{r ll t)8(r 2 ,t)---9{r Nl t)) } , (3) 

with (•••)/ denoting an average with respect to realiza- 
tions of the stationary forcing f(r,t) and of the velocity 
field u(r,t). It had been known for a long time that for 
high Re and Pe these functions are homogeneous func- 
tions of their arguments, i.e. 

F^\\r 1 ,...Xr N ) = \^F^\r 1 ,...r N ) , (4) 

with C^r being a scaling exponent that in general is 
anomalous, i.e. cannot be guessed from dimensional con- 
siderations. But only recently it became clear how these 
exponents are determined by the dynamical processes. 



To understand the progress made, we rewrite the dy- 
namical equation for 9(r,t) in the shorthand form 



d0(r,t) 
dt 



= L9(r,t) + f(r,t) 



(5) 



where in the present case C = u ■ V — re A. In recent work 
0, H it was clarified why and how passive fields exhibit 
anomalous scaling, when the velocity field is a generic 
turbulent field. The key is to consider a decaying problem 
associated with Eq. JSJ, in which the forcing f(r,t) is 
put to zero. The problem becomes then a linear initial 
value problem, 



89/ dt 
with a formal solution 



C9 



with the operator 



R(r,r',t)=T+exp[[ dsC{s)} 
Jo 



(6) 



(7) 



(8) 



and T + being the time ordering operator. Define next 
the time dependent correlation functions of the decaying 
problem: 

C {N \ ri ,---,r N ,t) = (9(n,t)---d(r N ,t)) . (9) 

Here pointed brackets without subscript / refer to the 
decaying object in which averaging is taken with respect 
to realizations of the velocity field and initial conditions. 
As a result of Eq. Q the decaying correlation func- 
tions are evolved by a propagator V^ N \r_, p,t), (with 
r = r 1 ,r 2 ,- ■ -,r N and p = pi,p 2 ,- ■ -,Pn) ■ 



C [N) {r,t) = dpV (N) (r,p,t) C^ l \p,0) 



(10) 



In writing this equation we made explicit use of the fact 
that the initial distribution of the passive field 0(r,0) is 
statistically independent of the advecting velocity field. 
Thus the operator V^ N '(r, p,t) can be written explicitly 

V [N \r,p,t) ee (R{n, Pl ,t)R(r 2 , p 2 ,t) • • • R(r N , p N ,t)) . 

(11) 
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The key finding 0,0 is that the operator V^ N \r, p, t) 
possesses a left eigenfunction of eigenvalue 1, i. e. there 
exists (for each N) a time-independent function (r ) 
satisfying 



Z (N) (r) = J dpZ (N \p)V (N Hr,p,t) . 



(12) 



The functions Z^ are referred to as "Statistically Pre- 
served Structures" , being invariant to the dynamics, even 
though the operator is strongly time dependent and decay- 
ing. How to form, from these functions, infinitely many 
conserved variables in the decaying problem was shown 
in p|, and is discussed again in Sect. III. The functions 
Z^ N \r) are homogeneous functions of their arguments, 
with anomalous scaling exponents Cn- 



(13) 



where ". . ." stand for sub-leading scaling terms. Since 
Eq. I|12l) contains Z^ N ' (r) on both sides, the scaling ex- 
ponent Cn cannot be determined from dimensional con- 
siderations, and it can be anomalous. More importantly, 
it was shown that the correlation functions of the forced 
case, F^ N ^(r) Eq. 10, have exactly the same scaling ex- 
ponents as Z( N > (r) [5j • In the scaling sense 



FW(r)~Z<*0(r) 



(14) 



This is how anomalous scaling in passive fields is under- 
stood. 

Besides exactly soluble examples in which the advect- 
ing velocity field is non-generic (i.e. ^-correlated in time) 
the existence of eigenfunctions of eigenvalue 1 of the time 
dependent propagator was demonstrated fully only in 
shell models of turbulence. While the present authors 
believe that shell models contain a lot of the robust prop- 
erties of real turbulence, this belief is not universally ac- 
cepted in the community. It is therefore necessary to 
demonstrate that the mechanism sketched above exists 
indeed in the full problem Eqs. Q-©- This had been 
done for 3'rd order correlations within the Lagrangian 
formulation in 0. The aim of this paper is to demon- 
strate this in the Eulerian frame, and for correlation func- 
tions of order 2, 4 and 6. 

In Section II we describe the simulations of Eqs. (JTJ 
and J2J an d the type of measurements that we performed. 
Sect. Ill is a theoretical digression, in which we analyze 
an exactly soluble simple model to guide ourselves as to 
how to analyze the numerical results to find the scaling 
forms of the nth oder propagators and their eigenfunc- 
tions. Sect. IV describes the analysis of the data, and 
Sect. V offers a discussion and a summary of the paper. 



II. SIMULATIONS 

We performed a Direct Numerical Simulation (DNS) of 
Eqs. and @ on a 2048 x 2048 2-D array. The forcing 



/ in Eq. Q is random, 5 correlated in time, isotropic 
and homogeneous. Its k depdendence is 

(|/(fc)|) oc k exp[-0.5(/c/1024) 2 ] . (15) 

This forcing is biased towards the small scales; this is 
done because of the inverse energy cascade that charac- 
terizes 2-dimensional turbulence. The fluid dissipation is 
modeled by a hyper-viscosity term proportional to A s u. 
In addition we employed a "friction" term proportional 
to u in order to stabilize the velocity field on the largest 
scales. The passive field 9 dissipates normally as shown 
in Eq. 

The simulations were performed for a decaying passive 
field 9, that is, the forcing / in Eq. J2J was put to zero. 
The initial conditions for the 9 field were of the form 



9(k, t = 0) = S(k- k )exp(i~f(k)) , 



(16) 



where 9(k, t) is the Fourier transform of the real space 
variable defined as: 



Z7T 



drh(r)e 



(17) 



(from now on we will omit the tilde above the functions 
and denote the k space functions only by their variables). 
7(fe) is a random variable in the interval 7 £ [0, 2ir] where 
7 (fe) = -7(-fe), insuring that 0*(-fc,O) = 9(k,0), and 
therefore that 6(r, t) is real. 

As the initial conditions for 9 and the forcing of the 
u fields are both homogeneous and isotropic so are the 
scalar correlation functions C^ N ^{r\ . . . r/v, t) defined in 
Eq. ©, at all times t. We measured the k space corre- 
lation functions: 



C^(k 1 ,...,k N ,t)S(k 1 
= (9(k u t)...9(k N ,t)) , 



+ k N ) 



(18) 



where again the average (...) is over initial conditions 
and over realizations of the u field. The delta function 
appears due to translational invariance. 

In accordance to Eq. (|T0|> we define the k space prop- 
agator: 

CW(fc,() = fdk!V^ N \k,k!,t) CW(fe',0) . (19) 



A. The 2-point propagator 

The only two point correlator that is not zero is 
C^(k,~k,t) = (6{k)9*(k)). Because of the isotropy of 
the initial conditions and the driving field, the correlator 
depends only on the magnitude of k. We can therefore 
consider the second order structure function S^ 2 '(k,t) — 
C (2) (fc, -k,t), and it's propagator P^(k,k',t). 

In discrete fc-space (as in our simulation on a grid) the 
propagator V^(k, kl , t) has a matrix representation. For 
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FIG. 1: P (2) (k, 300, t) for 10 different times. The first time is 
T30, i.e. the eddy turn over time at scale k = 30, in lightest 
grey. Later times, in darker and darker greys, are at 2t3o, 
3t3o, etc. until 10r 30 . 

the choice of initial conditions as in Eq. eg, SW{k,t) 
is simply the fcoth column of the propagator : 

5 (2) (k, t) = J dk'V (2) (k, k', t) S {2) (k', 0) = V {2) (fc, k , t) . 

In Figure ^ we plot such a column of the 2-point propa- 
gator, at 10 different times. Two properties of the prop- 
agator should be noticed: as time progresses the overall 
amplitude decreases due to the decay, while the maxi- 
mum moves from the initial ko to lower values of k. In 
Sects. Ill and IV we will find the scaling form of this 
propagator. 




10 1 10 2 
k 

FIG. 2: V (4) {k, 300, i) (upper panel) and T (6) (k, 300, t) (lower 
panel) for the same 10 different times as in Fig. 



B. Multi-point propagators 

In the case of the multi-point correlator the overall S 
function and the isotropy of the fields do not reduce the 
dependence to a one variable function. The propagators 
are therefore functions of many &;-vectors. It turns out 
however that measurement of the statistics for objects 
depending on many fc-vectors are very taxing. We opted 
therefore the extract from the DNS partial information 
on the dynamics of the 2N order structure functions with 
N > 1: 

S^ N Hk,t) = (\9(k,t)\ 2N ) . (21) 
Accordingly we define the reduced propagators : 

S^ 2N Hk,t) = Jdk'V^ N \k,k',t) S^ N Hk',Q) .(22) 

In Figure [2] we present the numerical results for of the 
fourth and sixth order reduced propagators, at the same 
ten different times. These objects display similar qual- 
itative behavior to that of the 2-point propagator. We 



will learn in the next section how to think about the 
scaling properties of these objects and how to re-plot the 
numerical results in proper rescaled variables. 



III. AN EXACTLY SOLVABLE CASE 

In order to motivate the data analysis presented in 
Sect. IV we turn now to the Kraichnan model for pas- 
sive scalar advection. The model is exactly solvable, and 
examining the analytic forms of the propagators offers 
clues to what should be expected in the generic case. In 
principle we could perform this analysis in terms of the 
Kraichnan model for a 2-dimensional passive scalar. In 
fact, it is sufficient to consider the shell-model version; 
the latter is very transparent to analytic manipulations, 
and for our purposes the results throw equally useful light 
on how to perform the analysis of the DNS results. 
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A. The Kraichnan shell model 

The Kraichnan shell model for passive scalar advec- 
tion as all other shell models of turbulent flows, is 

written in terms of Fourier components of the field. The 
Fourier components are restricted to shells denoted by 
the index m, and the equation takes on the form 



dt 



= C 



rn.n^n 



(23) 



£"m,n — ^m+l^m+l^m+l^ ^^rn^ m ^m— l,n ^k rn 8 7nn ■ 

Here k m are the shell k vector, k m = k X m for some ko 
and A. The shell components of the velocity field, u m (t), 
are Gaussian random variables, ^-correlated in time, sat- 
isfying: 

{u n (t)u* m (t')} = c 6 n , m d(t - t') A-« m . (24) 
Here £ is the scaling exponent of the u field. 

B. The 2nd order propagator 

The second order correlator satisfies the explicit equa- 
tion (see 0): 



dt 



Sg\t) = (Mi%-Kk 2 n )sW(t), (25) 



(2) 

Sn (t) is the shell equivalent of the 2nd order structure 

(26) 



function in Eq. (|2Ti|) . 

S^(t) = (\9 n (t)\ a ) . 
The time evolution operator has the form 

M nli = («n + OL n+1 )8 n , m 



(27) 



a n S r , 



On+iS. 



n+l,m 



where 



-c k 



<2 



-c k 2 X^ n , 



(28) 



and £2 = 2 — £ is the dimensional scaling of the 2-point 
correlation function. The operator M^ 1 has the follow- 
ing useful scaling property: 



\ -C2P nf( 2 ) _ M (2) 



(29) 



The second order propagator has, in the limit of vanish- 
ing viscosity, the explicit form: 



(2)> 



(30) 



Note that the time evolution operator and the propa- 
gator are both Hermitian, and thus admit an eigenvec- 
tor decomposition. The time evolution operator has two 



types of eigen-vectors ipn' q ^ and ipn ,q ^ which we can re- 
gard as slow modes and fast modes respectively. Here 
the index k stands for the eigen-mode index. The fast 
modes are dominated by the viscous term, their support 
is in the viscous range, and they can essentially be taken 
to be ipn'^ = <5„. 9 for some shell q with the shell in- 
dex q > rrid, above which viscosity dominates Eq. I|25|l . 
The transition shell is determined by the condition 
nk 2 — cok& , . For the fast modes we have the following 

• •I'd. ''I'd 

approximate Equation: 



2 Ji( 2 ><?) 



dt 



-nk n ipn 



(31) 



Kikq . 



therefore their eigen- values are 8 q 

For slow modes, ip^ ,q \ which have their support on 
shells smaller than m<j, the dissipative term can be ne- 
glected and they satisfy an eigen- function equation of the 
form: 



8 ib {2 ' q) - M (2) ^ 2 ' q) 

l J qYn — lvl n,mVm 



(32) 



For a sufficiently large inertial range we can use the 
scaling property 1(25) 1 to get: 

/Vi 2 ' 9) = \~ C2P M [ X. m+P ^ q) ■ (33) 
Shifting indices we rewrite, 

Mfti = X-^M^l] . (34) 
Define now a vector 

Xn = Vi-J , (35) 

we have: 



A HqXn — lvl n ,mX-n 



(36) 



We can thus define ^n' q+P ^ = Xn, an eigen-function 
of the time evolution operator with eigen-value 

X C2P 8 q . 

We can thefore conclude that the eigen-vectors and 
eigen-values may be obtained from each other by shift of 
the indices: 



fn+p 

Pq+P 



^(2,9) 



(37) 
(38) 



Fig. O demonstrates that indeed two different eigen- 
functions coincide once shifted with respect to each other. 
It follows from J3E1) that 



(39) 



In Fig. 0] the spectrum of the eigen-values is plotted, 
showing the two expected regions: a region for which the 
scaling of the eigen-values is oc k^ 2 for the slow modes, 
and tx fc 2 for the fast modes. The lowest eigenvalue 
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FIG. 3: A plot of the eigen- vectors ipn' 15 (squares), and ip 
(triangles), for n £ [0,30] 
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FIG. 4: The eigen-values, plotted from small to large. The 
eigen-values are well fitted by the analytic predictions 
(dashed line) for the slow modes, and by fc„ (dotted line) for 
the fast modes. The transition to viscous range occur at shell 
m d = 80 



is proportional to oto. The eigen-function ip n 2 '°^ associ- 
ated with it can be calculated explicitly, and it exhibits 
a normal, dimensional, scaling in the inertial range: 



^ 2 -°) oc (k n /k )-^ 



(40) 



Using Eqs ,l|57)l and gUJ, that imply that ipn'^ has a 
scailing 'tail' with exponent £2, starting after shell q we 
establish: 

ij) n 2 ' q) oc (fc„/fc g ) _C2 oc a q /a n q < n < m d , (41) 

We now use these results to learn how to re-scale our 
numerical data. Suppose that we started with an initial 

(2) 

condition in some shell n in the inertial range, 
S PlV ,. The fast modes do not contribute significantly to 
the time dependence, since they decay fast and anyhow 
have no support in the inertial range. Thus, using only 



the slow modes and Eqs. (|27|) . i|30|) > we have for p < n, 



q=l 



9=1 



V {2) (t) 



(42) 



In going from the first to the second line we inserted 

the 5-function initial conditions, getting thus a column 

( 2 ) 

of the propagator V^lt). Remember that in our DNS 
calculations we fixed the column index, and are interested 
in the scaling behavior with respect to the row index 
and time. Therefore we compute now (up to an overall 
dimensional constant), 



q=l 



p-m 

— V a 2 e a «+ mt = — 



f q— m+1 



pw( t )_n_p( 2 ) ( t ) 

p\n V / m\n \ f 



(43) 



where we have used the explicit form of the propagator 
in Eq. (|42|l and Eq. For m much smaller than p we 

can neglect the second term. Then we conclude that the 
propagator is a homogeneous function of the variables k p 
and t. Explicitly, multiplying by t, 

tV (2 \k p ,k n ,t) = (X-^ m t)V^(X m k p ,k n ,\-^ m t) (44) 

For fixed k n this is a homogeneous function of two argu- 
ments which can be always written in the form 



a n t 



(45) 



for some function A(x). The symmetry between p and n 
is restored by realizing that the asymptotic form of A(x) 
is 1/x for large x. 

This scaling form was also found in the study of the 
generic shell models of passive scalar advection (arbitrary 
time dependence in the velocity correlations) |5| and we 
will show that the same scaling form survives when we 
go back to the generic model studied by DNS. 



C. The multi point propagators 

The invariance of Eq. I|24|) under a uniform phase 
change n — > e^0 n dictates that the only non-zero cor- 
relation functions will have an equal number of variables 
n and conjugated variables 0* m . We define: 

c[™l im , h ,..., jN {t) = (0 il (t)...0 iN (t)0* 1 (t)...0* N (t)) 



(46) 
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As in the 2-point case we can define the respective multi- 
point differential time derivative operator, analogous to 
Eq. Ij25(l . In the limit of vanishing viscosity: 



dt -'-L hiji * >j v 1 

The respective propagator is defined by: 

^ N )a-\ - -n( 2Ar ) (+\n( 2N ) 



(27V) 



id\iL,j' 



1,3)1 ,3_ v 1 1 i 3_ v ' 

exp(tAf< 2JV >) 



i>ili')i' 



(47) 

(48) 
(49) 



Because of the time reversibility of the statistics of the 
u fields, and the anti-hermiticity of the un-averaged dif- 
ferential operator C in Eq. (|23|) : C* n n = —C niTn , both 
the time derivative operator M^ 2N ^ and the propagator 
itself are hermitian, and therefore admit eigenvalue de- 
composition. Furthermore the operator M^ 2N ^ has the 
following scaling property: 



M 



(2N) 



l+P-3 +P \i_+P, 3 +P 



(50) 



As in the two point case the dynamics within the scal- 
ing range is determined by slow modes. 



(2N,l) , (2N,l,k) 



h3 



M (2N) ,(2N,l,k) 



(51) 



In this equation the index I stands for a family of eigen- 
modes, and k for their index within the family. Each 
I family can be obtained from any one of its members 
by shifting the indices. The eigen-values of the modes 
within a given family can be obtained also by a shift : 



,(2N,l,k+p) 

^i+p,j+p 



(2N,l,k) 



R {2N,l) _ xC2PR (2N,l) 
Pk+p — A Pk 



(52) 
(53) 



Note that in the 2-point case we had only one family of 
eigen-modes. Here we added the index I to the eigen- 
values and eigen-modes to distinguish the different fam- 
ilies. Note also that in Eq. (j53(l the scaling exponent of 
the eigenvalues is (2; this stems from the scaling prop- 
erties of the differential operator M^ 2N \ cf. eq. (|50|l . 
However, the eigen-modes display in general anomalous 
scaling which can be represented by 



(2N,l,k) 



A 



-p (2 



(2N,l,k) 

i+pj+p 



(54) 



where &n.i is the (anomalous) scaling exponent of the 
Zth family of eigen-modes. 

Since the DNS were analyzed in terms of fused objects, 
we focus on initial conditions in Eq. (|4H|> for which all the 
2N indices are the same. We also measure the resulting 
structure functions 



q(2N) 



\2N\ 



(55) 



This procedure will extract a fused propagator for which 
there are only two indices, and we denote it bellow as 



^p\n W ' This propagator is not Hermitian, it has in 



general no eigen-function decomposition, but we can un- 
derstand the scaling form of any of its columns from the 
knowledge of the full propagator and its eigen- functions. 

The equivalent of Eq. g5J> for Sf N) {t = 0) = 5 k , n , 
p < n, summing over all families of slow modes (i.e. over 
the index I) is: 



I q=l 



E 



C, 



P 2<2 



T J2 a ^ e^ = P^>) .(56) 





l (a p a n ) 5 g=i 

In going from first to second line we have used the fact 
that 



,2N.l,p 

i>i.i... .i ^ 



-C,2N,l/(,2 



(57) 



This follows directly from Eqs. and 

Using the same argumentation as in Eq. (|43|) we get: 



V 



(2N) 



(*) 



E 



c, 



(a n t) & 



■A^ N ' l \a p t), (58) 



here again the C/'s are dimensional constants and the 
functions A( 2JV, ^(x) are some functions with the asymp- 
totic form A( 2N ' l \x) w x ~ for x large. We note that 
for sufficiently long times and if the scaling exponents are 
well separated, we can expect the sum to be dominated 
by the leading scaling exponent. 



IV. DATA ANALYSIS 

The detailed analysis that was possible for the simple 
model of Sect. Ill is not available for the generic model 
Eqs.QJ, (J2J. Our strategy is to assume that the scaling 
forms derived in the last section are still valid for the 
generic model, and to demonstrate this by replotting the 
data accordingly. We will see that these predictions are 
born out by the data. 



A. analysis of the 2-point porpagator 

In light of Eq. (|4*4*|) . we expect the 2nd order propa- 
gator fi( 2 \k,ko,t) to be, for a fixed fc , a homogeneous 
function of the variable kt 1 ^ 2 , and to decay as 1/t: 

V {2) {k,k ,t) oc -H^f^kt 1 ^ 2 ) , (59) 

where is some function. We test the correctness of 
this form in Fig. [S] To this aim we re-plot that data 




10" 10 



FIG. 5: tV (2) (k, 300, t) for 10 different times (light grey ear- 
liest dark grey last) 



shown in Fig. ^in different coordinates, multiplied by t 
and as a function of 



kt 1 '^ 



(60) 



The quality of the data collapse appears to strongly sup- 
port the proposed scaling form. We note that the data 
collapse is superior on the right of the maximum, and less 
convincing at its left. We believe that this stems from two 
reasons. First, there is better statistics for the right part 
of the curve, simply because it belongs to larger k vectors 
where the angular average is more extensive. Second, the 
left part of the curve is more sensitive to the finite size 
effects, particularly to the fact that the driving velocity 
field loses its scaling form close to L. With all the limita- 
tions of the numerical simulations we consider the data 
collapse as very satisfactory. 

If the prediction for the 2-point function holds, then 
we will have for the time dependent integral: 



V {2) {k,k',t)dk 



1 



H {2) {k)dk oc 



1 



(61) 

In Fig. H3 We show that indeed after an initial period 
the integral settles on a scaling form consistent with Eq. 

Using the form of the propagator we can establish the 
existence of a left zero mode of the 2-point propagator. 
Integrating over the two sides of Eq. Q12[l we have a time 
independent expression on the left hand side. We there- 
fore expect that the weighed integral of the propagator 
with the function will be constant: 



Z^(k)V {2) (k,k',t)d 2 k 



H{k) 



(62) 



dk = const 



10" 
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FIG. 6: The integral of the 2-point propagator, 
/ °° V {2) (k, 300, t)d 2 k, (solid line) as a function of time. 

The dotted line is the expected t~^ 1+1 ' c - 2 \ (time is measured 
in this figure, and in the ones that follow in units of T30 the 
eddy turnover time of k = 30) 
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FIG. 7: The integral of the 2-point propagator P (2) (k, 300, t) 
weighed by the left zero mode Z^(k) (solid line), compared 
to the integral of the un-weighed second order structure func- 
tion in dotted line. 



where we have used the fact that in an isotropic 2- 
dimensional systen we have Z^ 2 \k) ~ F^- 2 \k) oc k~ 2 ~& . 
We note that the constancy of this integral should be 
judged on the background of the decaying function, as 
done in Fig. [7| We see that while the 2nd order struc- 
ture function decays over three orders of magnitude, the 
"constant" objects changes by a factor of 2. The lack 
of constancy can be attributed to the sensitivity to the 
outer scale as seen in the data collapse in Fig. JSJ. If 
the collapse on the left side of the curve were perfect, so 
would be the constancy of the weighted integral. Note 
that in the calculation we have employed (2 = 0.67 in 
agreement with 
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FIG. 8: t^P (4) (fc, 300, t) (upper panel) and t& P (6) (k, 300, t) 
(lower panel) for 10 different times (light grey earliest dark 
grey last) 
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FIG. 9: The integral of the 4-point (upper panel) 

and 6-point (lower panel) propagators in solid line, 

J °°P (2m) (fc,300,t)d 2 fc, compared with the expected 
i -(i+Cajv)/Ca j in tne dotted line. 



B. analysis of the multi-point propagators 

Examining Ea.l jftSjl . we expect that for long times the 
functional form of the propagator for the fused correla- 
tion function defined in Ea. l|22|l is: 



V {2N) {k,k\t) oc 



-H^Xkt 1 '^) , (63) 



where C^n is the leading scaling exponent for the 2A'th 
correlation function. In Fig. |8jwe demonstrate the data 
collapse obtained by assuming this form for the 4'th and 
6'th order propagators. We notice the same excellent 
collapse on the right hand side of the function as in the 
2nd order propagator. Also the problems with the outer 
scale show up in a similar manner, giving less than im- 
pressive collapse of the left hand part of the function. 
For the present data collapse we employed simple scaling 
C,2N = N(2', our data does not support strongly anoma- 
lous exponents. 

Using the form of the fused multi-point correlators we 



can again predict the time behavior of their integral: 

1 



r {2N) (k,k',t)dk 



£(1 + C2n)/C2 



H(k)dk 



1 



oc 



£(1+Catf)/C2 



(64) 



In Fig. E|the time dependence of the integrals is plotted, 
showing agreement of similar quality to Fig. As in the 
2-point case, weighing the fused multi-point functions by 
the appropriate multi point fused steady state correlators 
should yield a constant: 



Z^ N \k)V^ N \k,k',t)d 2 k 



,oo £f(2A0( W l/Ca) 

oc ; . „ , , > — i-^kdk = 



^ + C2NfC2N I C2 



H( 2N 'Hk) 

fc!+C2iv 



(65) 



dk = const , 



where Z 2N (k) oc k~ 2 ~& N . The quality of the constancy 
is demonstrated in Fig. ^| Again finite size effects lead 
to some decrease in time of these weighted objects, which 
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FIG. 10: The integral of the 4-point propagator P (4) (k, 300, t) 
(upper panel) and the 6-point propagator (k, 300, t) 
weighed by the left zero modes Z w (k) and Z (6) (k) respec- 
tively in solid line, compared to the integral of the un-weighed 
object in dotted line 

nevertheless is very much reduced compared to the de- 
caying correlation functions. 

V. CONCLUDING REMARKS 

We have demonstrated in this paper that the generic 
advection of a passive scalar by a velocity field that obeys 
the Navier-Stokes equations can be discussed in terms of 
Eulerian Statistically Preserved Structures. By initiat- 
ing a decay with delta-function initial conditions (con- 
centrated on k — 300) we have found numerically the 
corresponding columns of the time-dependent propaga- 
tors for the 2nd, 4th and 6th order correlation functions 
(where for the 4th and 6th order objects we considered 
partial ("fused") information). Note that in contrast 



to the Lagrangian formulation of Statistically Preserved 
Structures [3|, for which there is no Preserved Structure 
corresponding to the 2nd order correlation, in the Eu- 
lerian formulation such an object exists and had been 
analyzed explicitly. We have used a simple (non-generic) 
model of passive scalar advection to guess the analytic 
scaling form of the propagators in the generic problem. 
The test for the relevance of this guessed form is the 
data collapse shown in Fig. 5 and 8. The guessed time 
dependence appears to be in close correspondence with 
the data as shown in Figs. 6 and 9. 

The analytic forms of the propagators predict the ex- 
istence of eigen-modes with eigenvalue 1 (Statistically 
Preserved Structures). We believe that this is the first 
demonstration of Eulerian Statistically Preserved Struc- 
tures in a generic flow. The numerical evidence for the 
constancy of the latter is encouraging, if not fully conclu- 
sive, as seen in Figs. 7 and 10. We attributed the (rela- 
tively small) decrease in amplitude of the putative Statis- 
tically Preserved Structures to the less than perfect data 
collapse at the largest scales (smallest k- vectors) that are 
seen in Figs. 5 and 8. These in turn stem from the inter- 
vention of the outer scale in our scaling range, a bound- 
ary effect that we did not succeed to eliminate in our 
modest-size simulations. It would be interesting to see 
whether larger 2-d simulations could remove this finite- 
size effect to demonstrate conclusively the constancy of 
the Statistically Preserved Structures. 

In concluding we wish to point out an additional ben- 
efit to the present formulation. Usually in modeling tur- 
bulent advection it is customary to resort to dubious con- 
cepts such as 'turbulent diffusion' in order to write a dif- 
fusion equation for the correlation functions. The present 
approach indicates a much better procedure, i.e. to find, 
for a given turbulent field, the form of the propagator 
which can be then used to provide analytic predictions 
for any initial (t=0) correlation function. This procedure 
may be quite attractive for more complex hydrodynamic 
flows with unusual boundaries or coherent structures. We 
believe that when the turbulence is sufficiently well de- 
veloped scaling forms for the propagator will exist, and 
once found can be used to solve efficiently any statistical 
initial value problem. 
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